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SUMMARY 


The  use  of  computers  to  optimize  free  parameters  of  a  system  has  become 
relatively  widespread  in  many  areas  of  engineering.   Parameter  optimization 
codes  have  been  written  for  that  purpose,  and  make  it  possible  for  a  design 
engineer,  once  he  has  developed  the  mathematics  of  a  system,  to  optimize  its 
parameters  according  to  some  criteria.   But  of  equal  interest  to  the  design 
engineer  is  the  sensitivity  of  the  optimized  criteria  to  departure  of  the 
parameters  away  from  their  optimal  value.   The  purpose  of  this  thesis  is  to 
show  ways  in  which  a  parameter  optimization  code  may  be  augmented  to  yield 
such  sensitivity  information. 

A  Fletcher  and  Powell  version  of  Davidon's  variable  metric  optimization 
search  technique  was  employed  to  optimize  multi-parameter  functions.   Their 
method  is  useful  in  that  it  computes  the  inverse  Hessian  matrix  (or  matrix  of 
second  derivatives) ,  which  completely  describes  the  curvature  of  the  function 
at  the  optimum.   Equations  were  developed  so  that  the  sensitivity  could  be 
expressed  in  a  meaningful  output  format.   This  was  made  possible  through  the 
use  of  matrix  inversion  and  eigenvalue  analysis  subroutines  which  were  ob- 
tained from  the  scientific  subroutine  library  of  the  IBM  360  91  and  used  in 
conjunction  with  a  digital  computer  code  employing  the  Fletcher  and  Powell 
technique.   Equality  constrained  optimization  problems  were  also  considered 
by  employing  the  penalty  factor  method  proposed  by  Courant  and  used  by  Kelley, 
Equations  analogous  to  the  use  of  Lagrangian  Multipliers  were  used  to  deter- 
mine the  cost  of  the  equality  constraint. 

Example  problems  are  offered  showing  the  optimal  solutions,  sensitivity 
data  at  the  optimum,  and  interpretation  of  that  data.   The  well  known 
Rosenbrock  function  was  used  to  exhibit  the  accuracy  of  the  methods  employed. 
A  typical  engineering  problem  was  solved  involving  the  sensitivity  of  an 
optimal  nuclear  rocket  engine  used  to  inject  a  payload  onto  an  interplan- 
etary trajectory.   The  results  indicated  that  the  thermal  power  of  the  re- 
actor and  the  ratio  of  length  to  diameter  of  the  core  could  be  varied  con- 
siderably from  their  optimal  values  with  little  cost.   The  power  density 
however  was  relatively  fixed  for  optimal  operation. 
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I.   INTRODUCTION 

The  use  of  computers  to  optimize  free  parameters  of  a  system  has  be- 
come relatively  widespread  in  many  areas  of  engineering.   Parameter  opti- 
mization codes  have  been  written  for  that  purpose,  and  make  it  possible 
for  a  design  engineer,  once  he  has  developed  the  mathematics  of  a  system, 
to  optimize  its  parameters  according  to  some  criteria.   But  of  equal  in- 
terest to  the  design  engineer  is  the  sensitivity  of  the  optimized  criteria 
to  departure  of  the  parameters  away  from  their  optimal  value.   The  purpose 
of  this  thesis  is  to  show  ways  in  which  a  parameter  optimization  code  may 
be  augmented  to  yield  such  sensitivity  information.   The  interest  in  this 
is  evident. 

Since  the  engineer  is  working  with  temperatures,  pressures,  masses, 
etc. ,  the  computed  optimal  solution  cannot  be  implemented  exactly.   Physical 
parameters  are  subject  to  uncontrollable  variations  and  the  resulting  de- 
parture from  the  optimum  may  be  quite  significant.   For  space  flight 
applications,  maximum  payloads  might  be  of  primary  interest  for  one  mission 
while  a  minimum  fuel  consumption  the  criterion  for  the  next.   Design  require- 
ments for  a  Venus  flyby  will  certainly  differ  from  those  of  a  Mars  lander 
or  Jupiter  probe,  and  yet  many  of  the  systems  must  be  flexible  enough  to 
be  useful  for  all  three.   Thus,  the  optimized  solution  must  also  contain 
significant  information  about  departures  from  the  optimum.   This  thesis 
addresses  that  problem;  i.e.  sensitivity  analysis  at  the  optimum. 

Parameter  optimization  algorithms  are  many  in  number  and  varied  in 
application.   They  differ  primarily  in  their  rate  of  convergence  and  the 
restrictions  imposed  on  the  function  under  consideration.   Nearly  all 
require  a  large  number  of  iterations  to  achieve  a  given  accuracy  in  locating 


the  optimum,  and  some  procedures  may  not  converge  from  a  poor  starting  point, 

Because,  near  the  optimum,  the  second  order  terms  in  the  Taylor  series 
expansion  dominate,  the  only  method  which  will  converge  quickly  for  a  gen- 
eral function  are  those  which  will  guarantee  to  find  the  optimum  of  a 
general  quadratic  function  speedily.   Fletcher  and  Powell  [1]*A  have  pro- 
duced such  a  method  which  was  based  upon  a  procedure  described  by  Davidon 
[2].   The  method  is  superior  to  other  techniques  which  possess  quadratic 
convergence  in  that  it  makes  use  of  information  determined  by  previous 
iterations  and  also  in  that  each  iteration  is  quick  and  simple  to  carry 
out.   The  primary  justification  for  its  use  however,  lies  in  the  fact  that 
the  method  yields  the  necessary  information  to  determine  the  curvature  of 
the  objective  function  at  the  optimum. 

The  method  of  Fletcher  and  Powell  estimates  and  continually  updates 
the  inverse  of  the  Hessian  matrix,  H,  (or  matrix  of  second  derivatives) 
during  the  optimization  search  so  that  a  close  approximation  to  the  true 
value  at  the  optimum  is  reached.  Through  the  use  of  the  eigenvalues  and 
eigenvectors  of  the  Hessian  matrix,  analysis  of  the  characteristics  of  the 
objective  function  in  the  neighborhood  of  the  optimum  is  possible. 

As  previously  mentioned,  it  is  of  interest  to  know  not  only  where  the 
optimum  is  located,  but  by  how  much  each  parameter  X.   may  be  changed  fr 


om 


*  Figures  in  parenthesis  indicate  references  listed  following  the  text, 
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its  optimal  value  before  a  significant  change,   A,   in  the  objective  func- 
tion occurs.   The  variation  of  the  x  .   may  be  independent  of  the  other  x  ., 
or  several  parameters  may  be  changed  simultaneously  in  a  co-ordinated  fashion 
away  from  the  optimum.   The  latter  variation  might  allow  for  significantly 
large  departures  from  the  optimum  for  a  specified  A  when  the  function 
possesses  the  characteristics  of  an  N-dimensional  "ridge". 

In  this  respect,  it  is  intended  that  this  thesis  may  be  used  to  evalu- 
ate and  analyze  the  optimum  of  any  general  function  of  N  independent 
variables  in  such  a  way  that  a  complete  sensitivity  at  the  optimum  is 
clearly  presented  in  a  useful  output  format.   It  is  shown  to  the  user 
which  of  the  specified  variables  may  be  changed  and  the  magnitudes  of  those 
changes  before  specified  degration  in  the  objective  function  will  be 
exceeded. 
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II.   PROBLEM  DEFINITION 

An  analysis  of  the  behavior  at.  the  optimum  of  an  N-variable  function 

is  possible  where  the  second  derivatives  are  known.   Suppose  that  Y   is 

a  real  valued  function  of  N  variables  with  continuous  first  and  second 

partials  and  possesses  a  relative  minimum  at  X  ,   then  the  first  deriva- 

o 

tive  will  vanish  and  by  Taylor  series  expansion: 


Y  -  Y    =AY=-AX 
opt        2   o 


8x2 


AX  +  Higher  Order  Terms 
o 


(1) 


where 


is  the  Hessian  matrix   (H)   of  Y  at  the  optimum, 


_°x  _ 


From  the  Taylor  series  approximation  (1)  we  find  that  the  gradient 


AY(X)  = 


JL 


(x) 


>x 


-1 


VY(X) 


(2) 


and  solving  for  X   yields : 


X  =  X  - 
o 


V2 

1  y 


n  -1 


(X) 


ax- 


VY(X) 


(3) 


so  that  if 


3x2 


1 


(X) 


were  known,  the  step  to  the  optimum  would  be 


given  by  (3) .   Some  optimization  algorithms  collect  information  which 


generate  an  approximation  to  the  inverse  Hessian  matrix  during  the  search 
for  the  optimum.   These  algorithms  provide  H~   so  that  the  problem  of 

obtaining   11,   in  the  use  of  such  an  algorithm  is  reduced  to  the  relatively 

92v 

simple  additional  task  of  inverting  the  N  x  N  matrix   — •u. 

J3x2_ 

Historically  the  method  is  similar  to  Newton's  method  which  minimizes  a 
function  Y(x),  x  on  N-vector,  by  generating  a  sequence  of  points  (compare 
with  equation  (3)) 


x(k+i)  .x0O.a0O  ptOO]-.  TOoo 


(k) 
where  a     is  an  appropriately  chosen  scaler  constant.   Fletcher  and 

Powell's  version  is  in  fact  a  quasi-Newton  method  which  uses  an  initial 

(k)  -  i 
estimate  and  computational  history  to  generate  an  estimate  to   [Hv   ] 

rather  than  performing  the  computational  work  of  evaluating  and  inverting 

the  matrix  at  each  step. 
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III.   DAVIDON'S  METHOD 

Davidon  introduced  a  variable  metric  method  which  was  the  first 
optimization  technique  to  use  past  information  estimating  the  inverse 
Hessian  matrix.   Fletcher  and  Powell  have  improved  upon  the  method  by 
simplifying  the  iteration  scheme  and  modifying  the  criterion  for  con- 
vergence.  Their  method,  which  numerically  determines  the  minimum  of 
a  function  of  several  variables,  combines  the  best  features  of  the  con- 
ventional gradient  method  and  Newton's  method,  namely  the  sureness  of 
convergence  of  the  former  and  the  quadratic  terminal  convergence  of 
the  latter.   An  excellent  exposition  of  the  method,  including  conver- 
gence proofs,  has  been  given  by  Fletcher  and  Powell.   For  the  purpose 
of  this  thesis  however,  it  will  suffice  here  to  state  the  algorithm  and 
point  out  the  usefulness  of  its  main  features  as  was  done  by  Kelley  [3] 

Let  us  denote  the  free  parameters  as  X.   and  the  function  to  be 

minimized  as   Y(x).   It  is  assumed  that  Y  has  continuous  first  and 

o 
second  partial  derivatives  with  respect  to  X.   Any  starting  point  X 

o 
is  chosen  according  to  some  criteria*.   At  the  starting  point  X  ,   the 

gradient  vector  Y   as  well  as  Y   itself,  is  evaluated.   A  change  is 

o 
then  made  in  X   according  to : 


AX  =  -a  H"1  Y  (4) 


*  The  freedom  in  the  choice  of  the  starting  point  depends  on  Y .   If  Y 
is  well  behaved,  this  choice  is  free.   In  other  cases,  the  starting 
point  must  be  close  to  the  optimum  to  assure  convergence. 


where  H"   is  a  positive  definite,  symmetric  matrix,  defining  a  metric 
in  the  X-hyperspace.   Its  initial  selection  is  otherwise  arbitrary.   For 
general  purposes  the  unit  matrix  may  be  used,  but  if  the  parameters  differ 
by  orders  of  magnitude  it  is  convenient  to  either  scale  them  or  to  estimate 
the  accuracy  with  which  they  are  to  be  determined,   a  >  0   is  a  step  size 
parameter. 

In  Davidon's  method,  the  one-dimensional  minimum  of  F  vs .  a  is 

O  -1 

obtained  along  the  vector  originating  in  X   in  the  direction  H   Y  . 

A. 

At  the  new  X,   the  gradient  vector  Y   is  again  evaluated.   The  H" 
matrix  is  updated  according  to 

H"1  +  AH"1  =  H  _1  +  A  +  B  (5) 


where 


T 

AX  AX 
A  =  

AX  AY 
x 


and  generates  the  inverse  of  H  in  N  steps  for  the  quadratic  function 
and  where 


-  H  _1  Ay  ay  t  h  _1 

B  =   x   x 

£Y  t  h  _1  ay 

x        x 


is  intended  to  cancel  the  initial  assumption  for  H  1    [4  ].   The  procedure 

*        o 

is  begun  again  with  the  new  values  of  X,   Y  ,  H 
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It  is  shown  in  [1]  that  H  remains  positive  definite  and  that,   as 
X  approaches  the  minimizing  point,   H~  approaches  Y  __  ~   evaluated  at 

XX 

the  minimum.  For  quadratic  Y,  in  N-dimensional  space,  the  minimum  is 
obtained  in,  at  most,  N  steps  (within  round-off  error):  the  method  is 
quadratically  convergent. 

For  more  general  functions  having  the  smoothness  properties  assumed, 
a  Taylor  expansion  through  quadratic  terms  provides  a  good  representation 
of  the  function  in  some  neighborhood  of  the  minimum.   With  H~   converged, 
the  minimum  of  Y  vs.  a  then  will  occur  for  a  =  1, 


AX  =  -a  H"1  Y 
x 


will  approach  the  value  given  by  Newton's  method,  namely  -Y     Y  . 

XX      X 
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IV.   SENSITIVITY  OF  THE  PERFORMANCE  FUNCTION  AT  THE  OPTIMUM 

Assuming  that  the  optimum  and  H    at  the  optimum  have  been  deter- 
mined and  that  H    has  been  inverted,  then  all  the  necessary  information 
has  been  collected  for  a  detailed  analysis  of  the  sensitivity  of  the  per- 
formance function  to  changes  in  the  performance  parameters. 

Let  the  criterion  function  be 


Y(xx,  x2, 


•  V 


and  assume  for  simplicity  that  the  origin,  i.e. 


x.,  ■  x_  =  =  x  =0 

12  n 


is  an  analytic  optimum  of  Y.   Around  the  optimum 


1   T 

Y  =  -j  X'HX 


(6) 


where 


H  = 


3^ 


3x.  3x . 
i   3 


s  a 


ij 


is  the  Hessian  matrix  of  Y  at  0. 


Independent  Variation  (one  at  a  time) 

If  one  parameter  at  a  time  is  allowed  to  depart  from  the  optimum, 
the  function  Y  is 


2  ix  i 


where  a..   is  the  corresponding  diagonal  term  of   H.   For  an  assigned 
change  A  in  Y  we  find 


±  Ax.  -  ± 

i 


2AY 


a.  . 


(7) 


where  ±  Ax.   is  the  allowable  change  in  x.   for  the  previously  assigned 
acceptable  variation  in  Y. 


Simultaneous  Variation 

If  one  allows  several  parameters  to  be  changed  in  a  co-ordinated  fashion 
away  from  the  optimum,  departures  from  0   for  a  specified   AY  may  be 
significantly  larger  than  those  shown  in  the  one-parameter-at-a-time  case 
To  illustrate  this,  consider  the  case  in  a  2-dimensional  parametric  space 
where  the  function  y   presents  the  characteristics  of  a  ridge  as  illustrated 
in  Figure  1. 

The  Y  =1A  line  intersects  the  x,   and  x   axes  at  distances  which 
are  far  less  than  the  values  of  x   and  x~   at  the  end  of  the  Y  =1A 
contour.   (Five  times  less  in  the  case  of  Figure  1). 
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A  Skew  Ridge  Illustrating  Advantage  of  Simultaneous  Variation 


Figure  1 
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Eigenvalues  and  Eigenvectors  of  _H 

To  be  able  to  analyze  the  characteristics  of  the  Y  surface  in  the 
neighborhood  of  0,  it  is  convenient  to  make  use  of  the  eigenvalues  and 
eigenvectors  of   H  at  that  point. 

It  is  assumed  here  that  the  final  use  of  this  analysis  will  be  in 
computer  programs,  and  that  matrices  eigenvalue  and  eigenvector  search 
subroutines  are  either  available  or  easily  implementable.   To  that  respect, 
we  note  that  H  is  a  symmetric,  positive  definite  matrix,  and  that  know- 
ledge of  this  fact  simplifies  the  implementation  of  such  subroutines. 

Let  X.   and  V.   be  respectively  the  eigenvalues  and  eigenvectors 
of  H  [5],  i.e. 

det  (H  -  XI)  =  0 


H  V.  =  X.  V. 
i    li 


Since  H  is  real  and  symmetric  the  eigenvalues  are  real  and  the 
eigenvectors  are  orthogonal  and  may  be  expressed  in  orthonormal  form: 


—  T   — 

V.   •  V.  =  1  (normal) 

i     i 


-T   - 

V.   •  V.  =  0   i  ^  j  (orthogonal) 


If  the  N-parameters  are  allowed  to  be  changed  in  the  direction  of  the 
ith  eigenvector,   i.e.   AX  -  k.  V  ,   then  the  degradation  in  the  perfor- 


mance 


function  AY  is: 
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1  .     -  T 

H  k.   \ 
1 


AY  =  -  k.   V.      H  k.    V.  (8) 

2      1      1  xi  v    J 


but 


__  _  _T     - 

Ii  V.    =  A.    V.      and     V      •   V     =   1 
ill  ii 


so   that 


2  i   i 
and  for  an  allowable  variation  A  in  Y 


AY  =  \   k.2  A.  (9) 


/2AY 
If  <10> 

k.   represents  how  far  in  the  direction  of  the  normalized  i    eigenvector 

one  may  travel  on  the  response  surface  before  degrading  the  performance 

function  by  AY.   It  is  observed  that  k  is  a  maximum  for  A  .  , 

mm 


so  that  the 

direction 

of  . 

least 

sensitivity 

to 

changes  in 

the 

parameters,   x. 

K 2 j 

the  eigenvector 

is  given 
with  A  = 

by 

A  . 
~  mm 

In  other  words  the  length  of  the  AX  vector  =   /  Z   AX.2  is  maximized  for 


X 


a  given  '  AY  in  the  performance  function. 

Consider  a  hypothetical  2-variable  optimization  problem  where  the  eigen- 
values of  H  at  0  assume  the  values  of  1  and  10,  i.e. 
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A  =1;   A2  =  10 


and  the  associated  orthonormal  eigenvectors  are; 


V1  =  (1,0)    and   V"2  =  (0,1) 


respectively.   The  contour  of  the  objective  function  at  the  optimum  will 
then  assume  the  shape  of  Figure  2. 

Observe  that  for  a  given  Ay,   the  distance  from  the  optimum  in  the 
direction  of  V..   is  considerably  greater  than  in  the  direction  of  V^. 
In  fact,  it  follows  directly  from  (10)  that  the  distance  is/  _2   greater. 

It  is  evident  that  the  relative  length  of  any  existing  ridges  in  the 
response  surface  will  be  determined  by  the  square  root  of  the  ratio  of 
any  two  eigenvalues  and  that  the  magnitude  of  the  sensitivity  for  a 
given  A  will  be  determined  by  the  square  root  of  the  eigenvalues. 

Let  two  distinct  solutions  of  optimization  problems  assume  the  values 
Solution  1  (as  before)  Solution  2 

Xu-1   Al2  =  10  ^-0.1   ^22  =  1 

Vx  =  (1,0)  V2  =  (0,1)  V±   =  (1,0)  V2  =  (0,1) 

The  response  surfaces  will  then  assume  the  shapes  of  Figure  3a  and  Figure 
3b  respectively. 

Notice  that  the  ratio  of  the  eigenvalues  has  remained  constant  and 
that  the  shape  of  the  response  surface  is  unchanged.   The  magnitude  of 
the  eigenvalues  was  decreased  by  a  factor  of  10  and  if  drawn  to  scale  the 
response  surface  would  be  expanded  in  all  directions  by  a  factor  of  /l0  . 
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x1  =  i 


A2  =  10 


Effect  of  Eigenvalue  Ratio  on  Response  Surface 


Figure  2 


-15- 


' — &-+ 


Au  -  1.0 


x12  -  10. 


T—  xl 


Figure  3a 


X' 


xn  -  0.1 


A22  -  1.0 


Figure  3b 


Effect  of  the  Magnitude  of  Eigenvalues  on  the  Shape  of  the 

Response  Surface 


Figure  3 
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Skew  Ridges 

The  previous  examples  have  pointed  out  that  when  the  matrix  H  has 
one  eigenvalue  much  smaller  than  all  other  ones  a  ridge  will  exist  in  the 
response  surface.   A  skew  ridge  exists  when  one  eigenvalue  is  much  smaller 
than  all  other  ones,  and  the  corresponding  eigenvector  is  not  parallel  to 
one  of  the  axis  in  X-space. 

If  a  ridge  is  parallel  to  one 
axis  in  X-space,  it  can  always 
be  removed  by  a  change  in  scale 
along  that  axis.   It  is  there- 
fore representative  of  the  way 
scales  are  chosen  rather  than 
a  characteristic  of  the  func- 
tion to  be  optimized. 

The  situation  is  completely  different  if  the  ridge  is  at  an  angle  to 
the  axis  because  no  change  in  scaling  can  remove  the  ridge.   Such  a  ridge 
reflects  an  interaction  between  parameters  in  the  way  in  which  they  affect 
the  criterion  function. 

Only  when  the  eigenvalues  and  eigenvectors  of  the  Hessian  matrix  are 
known,  can  such  a  ridge  be  discovered;  and  only  then  can  the  characteristics 

of  the  optimum  be  determined.  The  complete  search  of  eigenvalues  and  eigen- 
vectors of  H  and  the  derivation  of  the  resulting  sensitivity  information 
is  contained  in  the  computer  program  described  in  the  following  section. 


-17- 


-18- 


V.   COMPUTER  PROGRAM 

The  program  described  in  this  thesis  makes  use  of  a  standard, 
Powell-Fletcher  type,  parameter-optimization  computer  software 
package  (6).  It  contains  a  multi-dimensional  optimization  algorithm 
similar  to  Davidon's  variable  metric  method  as  was  previously  described. 

The  matrix  inversion  and  eigenvalue  analysis  subroutines  were  obtained 
from  the  scientific  subroutine  library  of  the  IBM  360  model  91.   A  listing 
of  these  subroutines  is  presented  in  the  Appendix.   A  simplified  flow 
diagram  illustrating  the  coordinated  use  of  these  programs  is  shown  on 
the  following  page. 

Necessary  Input  for  Sensitivity  Output 

The  input  is  the  expression  relating  the  objective  function  Y  to 

the  N  independent  variables  x..   The  expression  may  take  the  form  of 

■     x 

a  single  equation,  e.g. 


Y  =  f  (x±) 


or  may  comprise  any  number  of  subroutines  as  long  as  the  objective  function 
Y  and  the  N  independent  variables  are  defined.   The  allowable  departure  A 
from  the  optimum  is  also  a  required  input  and  must  be  specified  by  the  user. 
It  may  be  expressed  as  a  percentage  change  in  Y,  e.g. 


Y  =  f(x) 


•     

POWELL  -  FLETCHER 
OPTIMIZATION  ALGORITHM 


X     Y      H"1 
OPT'   OPT' 


MATRIX  INVERSION 


H 


■  '■ 
EIGENVALUES, 

EIGENVECTORS 

H,  A.,  V. 

SENSITIVITY 
EQUATIONS 


V 


AX.     (independent) 

(AX) .   (along  eigenvector) 


Flow  Chart  of  Computer  Programs 
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Ay  =  1%  y 


opt 


or     Y  =  99%  Y 


opt 


or  as  a  fixed  quantity,  e.g, 


AY  =  10. 


Output  Format 

Given  the  required  input: 


Y  =  f(x.) 


Y  ■  A    (a  constant) 

the  sensitivity  information  is  presented  in  the  following  table  for  con- 
venience in  analysis. 


Y        AY 
opt 


TABLE  I 


x 


X- 


X, 


"opt 


X 


N 


Ax  for  Specified  Ay 
Independent    Simultaneous  Variation 


Variation 


±Ax„ 


±Ax, 


±Ax 


N 


min 
Ax, 


Ax, 


Ax. 


X   = 


X 


max 


Ax. 


Ax, 


Ax 


N 
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VI.  NUMERICAL  EXAMPLE:   ROSENBROCK  FUNCTION 

In  order  to  illustrate  the  usefulness  and  accuracy  of  sensitivity 
analysis  a  standard  test  function  is  offered  as  an  example.   The  function 


(Xl,  x2)  =  100(x2  -  x^r  +  (1  -  xx> 


was  proposed  by  Rosenbrock  and  is  interesting  in  that  it  possesses  a 
steep-sided  ridge  following  the  curve  x   =  x   as  shown  in  Figure  4. 

The  exact  solution  of  the  problem  is  offered  along  with  the  computed 
values  so  that  a  comparison  can  be  made  and  the  effectiveness  of  the 
techniques  employed  may  be  evaluated. 

Problem:  minimize  Y  =  100 (x  -  x  2)2  +  (1.0  -  x  )2 

Exact  Solution: 


Y    =  0.0 
opt 


x±       =  1.0 


x2   =  1.0 


H  = 


802     -400 


-400 


200 


H 


-l 


0.5     1.0 


1.0     2.005 


where  all  numbers  given  are  exact. 

Fletcher  and  Powell's  version  of  Davidon's  variable  metric  technique 
found  an  optimum  solution: 

Y    =  10 "6  -  0.0 
opt 

x    =  1.00007  -  1.0 
x.,   =  1.00017  ~  1.0 


x2 


-2.0 


Start 


-1.0  J- 


Rosenbrock's  Curved  Valley  Function 


Figure  4 
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and  the  approximated  inverse  Hessian  matrix 


H' 


0.5055     1.0051 


1.0051     2.0035 


-l 


which  represents  an  error  of  less  than  1%. 

The  matrix  inversion  subroutine  with  the  approximated  H  ~A  as  its 
input  resulted  in  the  Hessian  matrix: 


H  = 


828.3 


-415.5 


-415.5 


209.0 


in  error  of  about  4%.   This  error  is  more  than  acceptable  for  the  purposes 
of  sensitivity  analysis. 


TABLE  II 


Sensitivity  Data 

of  Rosenbrock 

.  function: 

Ax 

for  Specified  Ay 

Y 

AY 

X. 

Independent 

Simultaneous  Variation 

opt 

opt 

Variation 

X    .      =   .399 
mm 

X2   -  1037. 

10~6 

0.10 

1.00007 

±.0155 

Ax  =  +.319 

/ 

Ax  =  +.012 

1.00017 

±.0309 

Ax2  =  +.633 

Ax2  =  -.006 

Analysis: 

If  the  response  surface  were  to  be  constructed  from  the  sensitivity 
data  shown  with  no  knowledge  of  the  function  under  consideration  it  would 
resemble  Figure  5. 
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2.0 


1.5  — 


1.0" 


0.5— 


(1.318,  1.635) 


Opt  (1.0,  1.0) 


(0.682,  0.365) 


0.5 


1.0 


1.5         2.0 


Predicted  Shape  of  Rosenbrock  Function  at  the  Optimum 


Figure  5 
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A  comparison  with  Figure  4  shows  that  the  sensitivity  analysis  at 
the  optimum  gives  the  desired   results.   The  above  figure  shows  that  the 
function  is  extemely  sensitive  to  changes  along  the  x-^  axis  (Ax  =  ±.015 
for   Y  =  0.1)   as  well  as  to  changes  along  the  x  axis  (Ax  =  ±.031).   In 
contrast,  it  is  shown  that  a  skew  ridge  exists  along  the  direction  defined 
by  Ax  =  .318,   Ax~  =  .633  and  that  if  varied  simultaneously  x   may 
be  changed  by  30%  and  x~  by  60%  before  Y  changes  by  0.1. 


-25- 


-26- 


VII.   EQUALITY  CONSTRAINTS 

In  the  optimization  of  multi-variable  problems  it  is  often  the  case 
that  only  certain  combinations  of  the  variables  are  either  meaningful 
or  acceptable.   The  imposed  restriction  usually  assumes  the  form  of  an 
equality  (or  inequality)  constraint: 


G(x  )  =  0 

i 


The  analytical  solution  of  such  a  problem  can  be  found  by  the  method 
of  Lagrangian  Multipliers  [7] ,  which  seeks  values  of  the  parameters  for 
which  the  modified  objective  function 

F*  =  F  +  AG  (11) 

is  stationary,  i.e. 


F*  =  F  +  AG  =0  (12) 

XXX 


Solving  for  A  yields 


F 
A  =  -  -r-  (13) 

G 
x 


A  is  meaningful  in  that  it  represents  the  cost  of  the  constraint  G.   If 
G  were  relaxed  by  1  unit  then  F  could  vary  by  X. 


In  the  general  case  however,  numerical  search  methods  are  unlikely 
to  locate  all  types  of  stationary  points  for  the  modified  function  using 
Lagrangian  Multipliers.   A  more  feasible  approach  is  that  of  penalty 
factors.  [8] 

Penalty  Factors 

In  the  use  of  penalty  factors  a  modified  function  which  incorporates 
the  constraint  is  defined  as 


F*  =  F  +  PG2  (14) 


where  P  is  a  large  positive-valued  constant  (for  minimization) .   If  P 
is  chosen  large  enough  then  G  is  forced  close  to  zero  in  the  search  for 
the  optimum.   At  the  optimal  solution  G   is  equal  to  some  small  quantity 
£.   If   £  is  not  within  an  acceptable  distance  from  G   then  P   is  in- 
creased until  an  optimum  is  found  which  satisfies  G  within   £.   The 
cost  of  the  constraint  can  be  found  in  a  manner  identical  to  the  method 
of  Lagrangian  Multipliers. 

At  the  optimum  of  the  modified  function  of  (14) 


F  *  =  F  +  2P£  G  =  0  (15) 

XXX 


so  that 


F 
2P£  =  -  -^  =  X  (16) 

x 


is  the  cost  of  the  constaint  G. 
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Although  the  idea  of  a  penalty  function  seems  to  have  been  conceived 
some  years  ago  (Courant  [1943]),  there  has  been  very  little  computational 
experience  with  regard  to  its  applications.   By  examining  equation  (16) 
it  is  observed  that  as   £  approaches  zero,   P  becomes  infinitely  large. 
Large  values  of  P  however  effectively  produce  a  sharp  ridge  in  the  con- 
tours of   F*,   and  most  search  techniques  are  troubled  by  the  existence 
of  such  a  ridge.   The  choice  of   P   is  therefore  a  compromise  between 
large  values  for  small  violations  in  G,   and  smaller  values  to  eliminate 
troublesome  ridges  in  the  modified  objective  function. 

In  using  Davidon's  optimization  technique  however,  it  was  discovered 
that  even  if  a  ridge  presents  no  problem  in  finding  the  optimal  solution, 
P  may  not  be  chosen  arbitrarily  large.   For  this  situation,   £  becomes 
so  small  that  changes  in  the  parameters  of  order  e     produce  corresponding 
changes  in  F  which  are  less  than  the  criterion  for  convergence  in  the 
search  for  the  optimum.   If  £  is  to  be  meaningful  it  must  be  large  enough 
to  possess  a  unique  solution.   In  other  words  changes  in  £  must  be  large 
enough  to  affect  the  terminal  convergence  of  the  optimization  search.   It 
is  therefore  necessary  to  possess  some  insight  into  the  problem  before 
the  penalty  factor  method  can  be  employed. 

We  may  note  that  for  any  given  equality-constrained  optimization 

problem, .  X     will  possess  a  unique  value.   Analytically  A   is  found  to 

be  2P£  as  £  approaches  zero  and  P  approaches  infinity.   Let  us 

denote  e     as  the  maximum  allowable  violation  in  the  equality  constraint 
max 

and   £  .    as  the  minimum  £  which  will  produce  a  unique  optimum  within 
min 

the  limits  of  the  convergence  criterion.   £  must  now  satisfy 


£  .   <  £  <  £ 

min        max 


for  a  meaningful  solution. 
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Values  for  X     and  P   are  estimated  within  one  or  two  orders  of 
magnitude  such  that 


2P,  „.  e  -  A,   N. 
(est)       (est) 


If  the  resulting  optimum  possesses  an  allowable   £   then  the  solution  is 
found  with  the  cost  of  the  constraint 


A  =  2Pe  (17) 


and  the  constraint  violation: 


£  =  G    -  G  (18) 

opt 


If  £  >  £     or  £  <  £  .    then  P  must  be  increased  or  decreased  re- 
max  mm  ■ 

spectively  until  an  allowable  £  is  found. 

Example  Problem 

To  illustrate  the  use  of  the  penalty  factor  method  for  equality  con- 
strained optimization  problems  the  following  two  examples  are  offered  for 
comparison  and  analysis. 


Let    F 


1(x1,  x2)  =  F2(xr  x2)  =  (xx  -  3)2  +  (x2  -  3) 


and    Gi^xi»  xo^  =  x l  +  x2  ~  *  =  ° 


G  (x  ,  x2)  =  x1  x2  -  4  =  0 
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X, 


■i h 


Constrained  minimum 


Figure  6a 


X  =   1.0 


G2  =  0 


Figure  6b 


-F 


Interpretation  of  Cost  of  the  Constraint  — - 


Figure  6 
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The  response  surfaces  are  shown  in  Figure  6,  where  the  equality 
constrained  function  has  a  minimum  value  at  the  same  point  for  both 
problems,  i.e.   x1  =  x~  =  2.0.   Since  the  unrestricted  function  is  the 
same  and  the  solution  lies  at  the  same  point,  then  the  gradient  F   is 
the  same  for  both  problems.   The  constraint  and  its  gradient  G^   are 
different  however  and  consequently  the  cost  of  the  constraint  A  will 
have  two  distinct  values. 

Analytically  it  can  be  shown  that  -F  /G  =  A  -   2     for  the  linear 

X   X 

constraint  x  +  x  =  4   and  that  -F  /G  =  A  =  1  for  the  hyperbolic 
constraint  x  x  =  4.   The  cost  of  the  constraint  has  been  reduced  in 
the  second  case  because  F  is  less  sensitive  to  changes  in  G~  at  the 
optimum. 

Computer  Results: 

For  F*  =  F  +  P  G  2  and     F*  =  F  +  P  G* 

and  P  =  100  P   =  100 

F*  =  1.980  F*  =  1.994 

x  =  2.005  xx  =  2.0009 

x2  =  2.005  x2  =  2.0020 

e  =  g*  -  G  =  .010  £  =  .0057 

A  =  1.988  A  =  1.141 

In  the  second  case,  a  minimum  was  found  which  was  very  close  to  the 
constraint,  e.g.   e  =  .0057.   As  a  result  A  =  1.141,   a  difference  of 
.141  from  the  analytical  value.  This  error  may  be  explained  by  examining 
the  magnitude  of   E.   An  E  of  .0057  is  very  small  for  the  problem  under 
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consideration.   Changes  in  the  variables  of  order   £  will  not  affect  the 

optimal  solution,  so  that   £   is  actually  less  than  £,,.,„  as  previously 

MIN     r        J 

defined. 

A  second  optimization  was  performed  with  P   reduced  to  50  in  order 
to  increase   £.   (Remember  that  it  has  been  shown  that   2Pe   is  constant), 
The  results  were: 


Y* 

= 

1.989 

Xl 

= 

2.0065 

X2 

= 

1.999 

£ 

= 

.011 

A 

=: 

1.07 

It  is  observed  that  decreasing  P  resulted  in  a  more  meaningful  value 
for  £  and  consequently  a  closer  approximation  to  the  const  constraint,   A, 
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VIII.   SENSITIVITY  ANALYSIS  OF  A  NUCLEAR  ROCKET  ENGINE 

As  was  previously  mentioned,  the  primary  purpose  of  a  sensitivity 
analysis  is  to  aid  the  design  engineer  in  constructing  a  system  which 
will  operate  in  an  optimal  fashion,  despite  variations  in  the  controlling 
parameters.   To  that  effect,  an  illustrative  example  is  presented  whereby 
the  engine  design  parameters  of  a  nuclear  rocket  are  optimized  in  order 
to  achieve  a  maximum  payload  at  a  specified  hyperbolic  velocity. 

A  set  of  mathematical  models  of  the  elements  of  nuclear  rocket 
engines,  suitable  for  detailed  systems  analysis,  has  been  developed  [9] 
which  constitutes  the  basis  for  a  digital  computer  program  called 
NUROC/SAC  (Nuclear  Rocket  _Systems  Analysis  Code) .   The  Code  has  been  used 
to  describe  a  number  of  existing  engines  and  the  results  obtained  were 
found  to  be  accurate  [10] . 

ESCAPE  [11]  is  another  computer  code  which  calculates  geocentric 
(or  planetocentric)  tangential  thrust  escape  trajectories  and  which  may 
be  used  in  conjunction  with  NUROC/SAC. 

Input  Parameters 

From  a  design  viewpoint  the  most  important  input  parameters  to 

NUROC/SAC  are: 

Q  thermal  power  of  the  reactor,  watts 

D  diameter  of  the  reactor  core,  meters 

L  length  of  core,  meters 

L/D  ratio  of  core  length  to  diameter 

T  maximum  allowable  core  material  temperature,  °K 

cmax 


p°  nozzle  chamber  stagnation  pressure,  n/cm2 

NAR  nozzle  area  ratio:   A   .  /A 

exxt   throat 


The  output  format  of  NUROC/SAC  consists  of  a  summary  of  operating  variables 
including  the  input  design  parameters,  defining  entirely  the  characteristics 
of  a  specific  nuclear  rocket  engine.   The  most  important  of  these  and  the 
ones  which  will  be  used  as  inputs  into  ESCAPE  are: 

m  total  engine  propellant  mass  flow 

F  total  engine  thrust 

m  total  engine  mass 

£i 

Additional  inputs  to  ESCAPE  which  must  be  specified  are: 

m  initial  mass  in  earth  orbit 

o 

H  initial  orbital  altitude 

o 

V  _  final  hyperbolic  excess  velocity  specified 

K„  tankage  to  propellant  mass  ratio 

The  output  of  ESCAPE  may  assume  a  variety  of  formats  but  the  payload 
delivered  at  the  specified  hyperbolic  velocity  is  the  payoff  in  the  opti- 
mization problems  which  follow.   The  payload  is  defined  as  the  initial  mass 
minus  the  engine  mass,  the  tankage  mass,  and  the  propellant  mass  required 

to  reach  V,  _. 
hf 

A  previous  study  [12]  has  indicated  that  within  the  range  described  by 

technological  constraints  the  payload  performance  of  the  nuclear  rocket 

will  increase  monotonically  with  the  maximum  allowable  core  material 
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temperature  and  the  nozzle  area  ratio,  and  will  decrease  monotonically 
with  the  chamber  stagnation  pressure.   Thus  choice  of  these  parameters 
is  limited  to  technologically  realizable  values.   Typical  values  presently 
in  use  are: 

T  =     3000  °K 

cmax 

p°         -     300  n/cm2 
NAR        =     100. 

The  power,  power  density,  diameter,  and  length  which  parameterize  the 
core  geometry  may  be  varied  in  a  constrained  but  optimal  fashion  to  describe 
an  engine  which  will  inject  the  maximum  pay load  into  a  specified  inter- 
planetary trajectory  for  a  given  initial  mass  in  earth  orbit.   The  problems 
which  follow  will  consider  an  initial  mass  of  100,000  kilograms  in  an  earth 
orbit  of  500  kilometers  and  will  optimize  the  core  geometry  to  inject  the 
maximum  payload  into  a  trajectory  described  by  a  hyperbolic  excess  velocity 
of  10,000  m/sec. 

Optimization  of  Core  Geometry 

It  is  interesting  to  note  that  if  any  three  of  the  four  engine 
parameters,  power  (Q)  ,  power  density  (ip) ,  diameter  (D)  and  length  (L) , 
are  specified  then  the  fourth  is  automatically  determined.   Since  the 
length  and  diameter  describe  the  volume,  the  power  density  can  be  ex- 
pressed as  a  function  of  Q,  L,  D,  i.e., 

2 
ip  =  4Q/ttD  L 

Thus  the  optimization  of  core  geometry  is  reduced  to  a  3-dimensional 
problem. 
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Maximize  Payload:   Free  Parameters   Q,  \p ,   L/D 

The  following  example  maximizes  the  payload  delivered  by  a  nuclear 
rocket  with  free  parameters:   power  (Q) ,  power  density  (ty) ,  and  ratio 
of  length  to  diameter  (L/D). 

Starting  Point: 


Q   =  x  =  2000.  MWt 

\p       =  x  =  6.0  x  10 9  W/m3 


L/D  =  x  =  4.0 


TABLE  III 


Sensitivity  Data 

at  the  Optimum 

Ax  for  Specified  Ay 

Y   . 
opt 

AY  =  1% 

x      : 

opt 

independent 
Variation 

Simultaneo 
X   =  39.5 

us  Variatior 
X   =  1105. 

i 
X   -  6150 

31,630  Kg. 

316.  Kg 

Q=1832  MWt 

AQ=±333. 

68. 

-28. 

28. 

ip=8. 84x10  9 
W/m3 

Aip=±3.03g 
(xlO  ) 

3.88 

(xl0y 

-.062  q 
i     (xlO  ) 

-•0179N 

(xl0yf 

L 

/D=4.80 

AL/D=±.57 

.387 

.677 

.199 

■ 

Analysis : 

When  both  \p     and  L/D  are  free  parameters  the  optimized  values 
become  so  high  that  the  cost  of  the  uranium  necessary  to  make  such  a 
reactor  critical  would  become  prohibitively  expensive.   Also,  a  length 
to  diameter  ratio  of  approximately  5  and  a  power  density  of  9x10  W/m 
would  describe  a  highly  inefficient  reactor  due  to  excessive  neutron 
leakage  through  the  core  ends.   Other  important  design  factors  (such  as 
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shielding),  which  are  not  now  contained  in  NUROC/SAC,  indicate  that  the  length 
should  be  only  slightly  larger  than  the  diameter.   It  should  also  be  noted 
in  the  sensitivity  data  that  \p     is  quite  flexible   (Ai|»  =  ±3xl09)   and 
that  smaller  values  may  be  chosen  with  little  resulting  loss  in  payload. 

This  problem  should  serve  as  an  example  that  one  cannot  randomly 
optimize  variables  or  undertake  a  sensitivity  analysis  without  first 
acquiring  a  knowledge  of  the  system  under  consideration.   With  these 
facts  in  mind,  a  second  optimization  problem  is  solved  in  two-dimensions 
with  L/D  fixed  at  1.5. 

Maximize  Payload:  Free  Parameters  Q,  ip 
Starting  Point: 

Q  =  x  =  2000  MWt 

ip  =  x2  =   6.0xl09   W/m3 

TABLE   IV 
Sensitivity  Data  at   the  Optimum 


opt 


AY  =  1% 


X 


30,938Kg. 


309. 


op' 


Ax  for  Specified  Ay 
Independent    Simultaneous  Variation 


Q=1823  MWt 

i|;=4.97xl09 
W/m: 


Variation 


AQ  =  ±339. 

A^  -  +.160 
(xlO9) 


X   =  5151 


340. 

-.02 
(xlO9) 


X  =  23630 


9.6 
.159 


9 
(xlO  ) 


Analysis: 

The  restriction  on  L/D   (fixed  at  1.5)  results  in  a  loss  in  optimal 
payload  of  only  2  percent.   It  is  interesting  to  note  that  the  optimal 
power  and  allowable  deviation  are  almost  identical  to  the  three-variable 
problem.   The  power  density  however  has  a  considerably  lower  optimal 
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value  and  A^   is  much  smaller  when  L/D   is  fixed. 

The  sensitivity  data  reveals  that  any  power  between  1500  and  2150  MWt 
may  be  used  if  ty     is  held  constant  with  only  a  small  (1%)  resulting  loss 
in  payload.   In  contrast  the  power  density  must  remain  close  (within  4%) 
to  5.0xl09  W/M3.   The  fact  that  the  eigenvectors  are  located  close  to  the 
Q  and  ty     axes  implies  that  little  more  can  be  gained  by  varying   Q  and 
Tp     simultaneously. 

Maximize  Payload:   Free  Parameters   Q,   L/D 

With  the  discovery  that  \p     and  L/D   cannot  both  be  free  and  having 
fixed  L/D  at  1.5  it  would  now  be  advantageous  to  fix  ty     and  optimize 

Q 

on  Q  and  L/D.   A  value  of  3.0x10  W/M   is  chosen  for  the  power  density 
based  on  accepted  values  for  existing  nuclear  rocket  engines. 


Starting  Point: 


Q  -  x  -  2000  MWt 
L/D  =  x2  =  2.0 


TABLE  V 


Sensitivity  Data  at  the  Optimum 


"opt 


AY  =  1% 


X 


opt 


Ax  for  Specified  Ay 
Independent     Simultaneous  Variation 
Variation         X   =  935.  X   =   3867. 


30,675  Kg. 


306. 


Q  =  1542  MWt 
L/D  =  1.67 


AQ  =  ±397 
AL/D  =  ±.78 


AQ  =  114 
AL/D  =  0.79 


AQ  =  390 
AL/D  =  -.056 


Analysis: 

When  the  power  density  is  fixed  both  the  power  and  the  ratio  of  length 
to  diameter  have  large  acceptable  variations.   The  power  may  vary  from  1150 
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to  1950  with  L/D  and   ip  fixed,  and   L/D  may  vary  from  1.9  to  2.45  with 
Q  and  ip  fixed.   This  allows  for  considerable  flexibility  in  designing 
an  optimal  engine  as  long  as  the  power  density  remains  relatively  constant. 

The  least  sensitive  direction  of  change  is  associated  with  the  mini- 
mum eigenvalue.   The  associated  eigenvector  indicates  that   L/D  may  still 
vary  from  .9  to  2.45  and   Q  need  not  remain  fixed,  but  may  vary  by  ±114  MWt 
as  long  as  the  direction  of  change  from  the  optimal  value  coincides  with 
changes  in  L/D. 


Verification  of  Results 

In  order  to  verify  that  the  predicted  sensitivity  of  the  optimized 
variables  is  accurate,  the  allowable  deviations  were  substituted  into 
NUROC/SAC  and  ESCAPE  and  the  resulting  payload  was  computed.   A  comparison 
was  made  to  check  if  the  payload  remained  within  the  predicted  1%  of  the 
maximum.   The  optimal  values  were  first  rounded-off  to  Q  =  1550  ±400  MWt; 
L/D  =  1.67  ±.75. 

TABLE  VI 
Comparison: 

Maximum  payload  minus  1%  =  30,370  Kg. 


Independent  Variation 
I.   AL/D,  Q  fixed 

II.   AQ,   L/D  fixed 


L/D 


Simultaneous  Variation 


Payload  (Kg) 


.90 

30,271 

2.45 

30,520 

£ 

1150.  MWt 

30,295 

1950.  MWt 

30,420 

9l 

L/D 

1435 

.90 

30,135 

1665 

2.45 

30,445 
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The  table  of  comparisons  indicates  that  if   Q  and  L/D  are  increased 
either  independently  or  simultaneously  the  loss  in  payload  is  even  less 
than  the  predicted  1%.   If  decreased  the  deviations  become  slightly  larger 
than  1%  but  are  still  highly  accurate.   This  may  not  always  be  the  case. 

It  must  be  realized  that  the  sensitivity  analysis  is  accurate  only 
at  the  optimum.   For  functions  which  are  very  flat  (i.e.  insensitive  to 
change)  the  predicted  variations  may  be  quite  large,  as  indeed  they  are 
in  the  example  given.   When  the  variations  are  substituted,  the  function 
is  no  longer  close  to  the  optimum  and  results  may  vary  considerably  from 
those  predicted.   For  this  reason  it  is  prudent  to  verify  results  especially 
at  points  far  from  the  optimum. 

Another  reason  for  verifying  results  is  that  there  exist  small  in- 
herent errors  in  the  optimization  search,  matrix  inversion  and  eigen 
analysis  subroutines.   Care  should  also  be  taken  in  defining  the  conver- 
gence criterion  because  if . the  computer  is  forced  to  make  repeated  searches 
near  the  optimum  the  values  for  the  inverse  Hessian  matrix  will  be  destroyed. 
Any  resulting  sensitivity  analysis  will  be  meaningless. 

Application  of  Penalty  Factor 

In  the  previous  two-dimensional  optimization  problem  an  optimal  power 
of  1550  MWt  and  length  to  diameter  ratio  of  1.67  was  established  for  a 
fixed  power  density  of  3.0x10  M/W  .   As  was  previously  mentioned  for 
equality  constrained  problems  it  is  of  interest  to  know  the  cost  of  the 
restriction  on  power  density  X.      The  following  example  employed  the 
penalty  factor  method  to  optimize  a  three-variable  problem  with  the  constraint; 


G(x3)  -  x3  -  3.0xl09  =  0 
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Starting  Point: 


Q  =  x  =  2000  MWt 

L/D  =  x2  =  6.0 

ip  =  x   =  4.0x10°  W/M3 

Penalty  Factor  =  P  =  10,000 

Results: 

Payload  =  30,670  Kg. 

Q  =  1547  MWt. 

L/D  =1.7 

^  =  3. 020x10 9  W/M3 

e  =  g*  -  G  =  .020 

X  =  2P£  =  406. 

The  cost  of  keeping  the  power  density  constant  is  406  Kg.  of  payload 
Per  (ip  W/mJ). 
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IX.   CONCLUSIONS  AND  RECOMMENDATIONS 

A  sensitivity  analysis  at  the  optimum  has  been  performed  through 
the  use  of  existing  computer  codes.   By  simply  defining  an  objective 
function  and  any  number  of  independent  variables  an  optimized  solution 
is  found.   In  addition,  it  has  been  shown  that  the  sensitivity  of  each 
of  the  variables  at  the  optimum  may  be  given  in  a  useful  format.   With 
all  others  held  constant,  the  range  that  each  variable  may  assume 
before  a  specified  degradation  in  performance  occurs  is  given.   Also, 
the  length  and  direction  of  the  axis  of  all  ridges  in  the  objective 
function  are  listed  in  order  of  decreasing  magnitude.   With  this 
information  at  hand  a  complete  knowledge  of  the  sensitivity  of  all  the 
input  parameters  is  available  to  the  user  and  decisions  regarding 
optimal  parameter  choice  can  be  made,  taking  into  account  the 
flexibility  given  by  the  sensitivity  knowledge. 

The  methods  employed  were  shown  to  be  highly  accurate  when  dealing 
with  purely  mathematical  problems.   The  sensitivities  of  Rosenbrock's 
parabolic  valley  function  were  found  with  little  difficulty.   Engineering 
applications,  on  the  other  hand,  require  a  great  deal  of  insight  into 
the  problem  before  a  sensitivity  analysis  may  be  started.   Scaling  the 
variables  is  important  so  that  the  sensitivity  information  is  meaningful. 
The  nuclear  rocket  example  studied  the  engine  performance  related  to 
power,  power  density  and  ratio  of  length  to  diameter  of  the  nuclear  core. 

A  scaling  problem  existed  because  the  power  in  watts  and  power  density 

9 
in  watts/cubic  meter  were  of  the  order  of  magnitude  of  10  ,  while  the 

length/diameter  was  approximately  unity.   Usually  such  a  problem  can  be 


avoided  by  normalizing  each  variable  by  dividing  it  by  the  order  of 
magnitude.   This  would  be  fine  if  the  acceptable  deviations  occurred 
within  the  same  number  of  significant  figures.   However,  if  one  variable 
is  much  less  sensitive  to  change  than  all  the  others,  then  the  eigen- 
vector associated  with  the  minimum  eigenvalue  will  be  dominated  by  that 
component  and  very  little  knowledge  can  be  gained  about  the  other 
variables  in  that  direction.   It  would  be  convenient  to  modify  the 
computer  program  so  that  the  sensitivities  would  be  normalized  instead 
of  the  variables.   Of  course,  the  user  could  always  eliminate  the 
problem  by  normalizing  the  variables  the  first  time,  and  after  observing 
the  resulting  deviations,  normalize  the  sensitivities  and  optimize 
again.   The  optimized  solution  will  be  the  same  but  the  sensitivity 
data  will  be  more  accurate  and  meaningful. 

Another  problem  experienced  when  working  with  engineering  problems 
was  in  defining  a  convergence  criterion.   For  the  Rosenbrock  function  a 
change  in  the  variables  of  10"^  produced  significant  changes  in  the 
objective  function.  When  optimizing  the  nuclear  rocker,  however,  the 
maximum  payload  was  essentially  determined  when  the  normalized  variables 
were  accurate  in  the  second  decimal  place.   With  the  convergence  criterion 
set  at  10"  ,  the  program  made  over  a  hundred  more  iterations  before 
stopping  and  found  an  optimum  which  was  only  about  one  kilogram  more  in 
payload.  When  the  optimization  search  stays  close  to  the  optimum  for 
many  iterations,  the  inverse  Hessian  matrix  is  destroyed  and  any 
resulting  sensitivity  analysis  has  no  meaning.   A  modification  in  the 
program  is  required  so  that  when  the  optimum  has  essentially  been 
determined,  the  inverse  Hessian  matrix  can  be  stored  and  any  subsequent 
refinement  of  the  optimum  will  not  affect  the  sensitivity  data. 
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Equality  constrained  optimization  has  been  discussed  but  not  in 
detail.   The  penalty  factor  method  used  by  Kelley  was  implemented 
to  determine  an  approximation  to  what  is  usually  referred  to  as  the 
cost  of  the  constraint. 

It  would  also  be  of  interest  to  know  how  far  in  the  plane  of  the 
constraint  and  normal  to  it,  one  may  travel  on  the  response  surface 
before  reaching  a  specified  change  in  the  optimum.   This  thesis  has  not 
considered  such  a  sensitivity  analysis.   The  techniques  involved  are 
similar,  but  more  attention  is  needed  in  this  area. 

Inequality  constrained  optimization  problems  are  generally  no  more 
difficult  to  solve  than  equality  constrained  ones.  For  well  behaved 
functions,  the  solution  either  does  not  violate  the  constraint  or  lies 
on  the  boundary  and  may  be  treated  as  an  equality  constraint.   In 
either  case,  the  techniques  developed  in  this  thesis  can  be  applied. 
Care  should  be  taken  to  try;.a  variety  of  starting  points  so  that  if  the 
same  solution  is  reached,  the  function  can  be  assumed  to  be  well 
behaved . 
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APPENDIX  A 

This  appendix  describes  the  implementation  of  the  variable  metric  op- 
timization method  and  matrix  inversion  and  eigenvalue  analysis  computer 
programs  for  the  purpose  of  sensitivity  analysis.   The  input  variables 
and  their  definitions  are  shown  along  with  a  flow  diagram  illustrating 
how  the  programs  were  modified  to  run  successively  and  produce  the  desired 
sensitivity  data.   A  typical  output  format  is  also  given. 

The  matrix  inversion  and  eigenvalue  analysis  subroutine  listings  con- 
tain brief  descriptions  of  the  methods  employed  and  their  references.   For 
a  complete  description  of  the  computer  coded  variable  metric  method  consult 
reference  6. 


BLOCK  DIAGRAM  OF  CONTROLLING  PROGRAM 


A- 2 


START 


LOAD 

N,  STEP,   a,   3, 
NEXP,  DELLAM 

i  CHANGE,  JPOINT 
X(J),  X(J)MIN, 

:  X(J)MAX,  DELX(J) 


TEST 
N  <  10 

X(J)  >  XMIN(J) 
X(J)  <  XMAX(J) 


FALSE 

HALT 


FALSE  jqjIN  ( J  )  =X  (  J  ) 
XMAX(J)=X(J) 


TRUE 


CALL  VARMET 

(N,  F,  Xs  XMIN,  XMAX, 
DELX,  STEP,  a,  6, 
NEXP,  DELLAM,  CHANGE, 

JPOINT,  JFLOW,  H) 


TRUE 


RETURN 


JFLOW 
N.E.2 


V 


FALSE 


PRINT  OPTIMUM  AND 
!  CO-ORDINATES 


:  FB,  X(J)  J-l,  N 


' 


TEST  CONSTRAINT 
G  =  0 


TRUE 


NGRD  =  0 
NIMP  =  0 


INITIAL  CALL 
MODEL  (N,F,X) 


/_ 


NEVAL 

FB 

JFLOW 


1 
F 

0 


INITIAL  PRINTOUT 
NIMP,  NGRD,  NEVAL, 
FB,  X(J),  J=l,  N 


.  PRINT  INVERSE 

;  HESSIAN  MATRIX 

H(I,J) 
J — 


CHANGE  STORAGE 
MODE 
HH(I)  =  H(I,J) 


V 


MATRIX 
INVERSION 


OPTIMIZATION  LOOP 


A- 3 


=  GS  -GG 


=  2P£ 


PRINT 


£,  X 


CALL  MINV 


(HH,  N,  D,  L,  M) 


RETURN 


CALL  EIGEN 


PRINT  HESSIAN 
MATRIX 
HH(I) 


CHANGE  STORAGE 

MODE 

H(I,J)  =  HH(I) 


STORE  DIAGONALS 
UNCOR(I)  =  H(I,I) 


j 
1 

■ 
MV  =  0 

i 

1 

>-\ 


(HH,  R,  N,  MV) 

RETURN 


- 

CHANGE 

STORAGE  MODE 

H(I,J) 

=  HH(I) 

STORE  ; 

EIGENVALUES 

DECEIG(I)  -  H(I,I) 

SENSITIVITY  EQUATIONS 


DY   = 


1  DX   =   /  1^1 

1    >/  ki 


2DY 


V.    =  k.  V. 
l       li 


EIGENVALUE 
ANALYSIS 


READY  FOR  OUTPUT 


PRINT 
OPTIMAL  Y   VECTOR  OF  OPTIMAL  X 
ASSIGNED  VARIATION  'DY' 

INDEPENDENT  VARIATION  ±DX. 

l 

SIMULTANEOUS  VARIATION  X.,  V.. 


EXIT 


Variable  Metric  Program 


Input  variables 

N 

STEP 

ALPHA 

BETA 

NEXP 

DELLAM 

CHANGE 
JPOINT 


X(J) 
XMIN(J) 
XMAX(J) 
DECX(J) 


and  definitions: 

the  number  of  independent  variables 

initial  step  size  for  X   components 

factor  by  which  step  size  is  increased  (a  >_  1) 

factor  by  which  step  size  is  decreased  (0  <_  fc  _<  1) 

limit  on  number  of  experiments  along  a  vector  line 

termination  criterion  for  vector  search  expressed  as  a 

range  fraction 

program  termination  criterion  expressed  as  a  range  fraction 

control  for  type  of  numerical  differentiation  desired  in 

gradient  subroutine. 

1  =   forward  difference 

2  =  backward  difference 

3  =  central  difference 

vector  of  initial  values  for  independent  variables 

lower  bound  for  each  independent  variable 

upper  bound  for  each  independent  variable 

step  size  for  each  independent  variable  used  in  gradient 

subroutine  for  numerical  differentiation. 


There  are  N  +  2  data  cards  required  as  inputs  for  each  program  execution. 
The  first  is  an  identity  card  written  in  any  format.   The  second  card  contains 
the  first  eight  input  variables  listed  above  in  the  following  format  fields. 
(110,  3F10.4,  110,  2F10.4,  110) 

Each  of  the  N  remaining  data  cards  contains  the  initial  value,  range 
and  step  size  of  the  independent  variables  according  to  the  following  format 
statement.   FORMAT  (4F15.5) 
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A  typical  output  from  z'r.ci   variable  metric  optimization  program  follows. 
It  should  be  noted  that  for  the  inverse  Hessian  matrix  to  be  printed  properly 
the  array  designated  H  must  be  dimensioned  exactly,  i.e.  if  N  -  2  then 
DIMENSION  H  (2,2). 
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)B0 


:ep 

,PHA 
IT  A 
!XP 
■LAM 

IANGE 
»OINT 


C.10C0Q 

1.00000 

C. 50000 

5 

0.01000 

0.0CC1C 

1 


1.20000 
1.0000C 


~  XMIN 

-5.00000 
-5.00000 


XMAX 

5. 00000 
5.COO00 


D2LX 

C.C010C 
0.00100 


NIMF 


NGRD       NEVAL 


FUNCTION 


INDEPENDENT    VARIABLES 


0.584002 
C.26200E 


01 

00 


-0.  120002 
;-C.51186B 


C. 680542-01 

0.10833E-01 

C.40146E-12 

1.384532-1 


or 

00 

-C.35741E  00 
0.142592  00 
0.59605E-07 
.S3192E-G6 
202302-06 


i  0.  1000CE  01 
-0.3576  32-06 
-C.22444Z  00 
0. 895462-01 
-0.417232-06 

.  11S362- 
-C. 315732-08 


TIMIZATICN    COMPLETE    APTEI-        "ST  FUNCTION    EVALUATIONS 

6  9  53  0.396642-13  -0.20233E-06 

THE    INVERSE    HESSIAN    MATRIX    IS 


-0. 325562-08 


0.4999 


HESSIAN    MATRIX 


3.9973 


*ERRCR***  EN- 

PRCGEAMME    WAS    EXECUTING    LINE 


65    IN    ROUTINE    M/PROG    WHEN    TERMINATION    OCC 
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Matrix  Inversion  and  Eigenvalue  Analysis 

The  example  output  of  the  variable  metric  computer  code  shown  previously 
lists  the  optimum  point  and  the  inverse  Hessian  matrix  at  the  optimum. 
Before  the  matrix  inversion  subroutine  (MINV)  may  be  implemented  the  upper 
triangular  elements  of  the  symmetric  H  matrix  must  be  stored  in  a  singly 
dimensioned  array  (HH) .   This  was  necessary  because  of  the  input  format 
utilized  in  the  calling  sequence  of  MINV  and  was  accomplished  by  inserting 
a  'DO'  loop  in  the  main  program. 

Once  the  matrix  has  been  inverted,  the  elements  of  the  Hessian  are 
printed  and  the  diagonal  elements  stored  for  calculation  of  sensitivity 
information.   Implementation  of  the  eigenvalue  analysis  subroutine 
(EIGEN)  merely  requires  defining  the  constant  MV,  i.e. 
MV  =  0  eigenvalues  only 
MV  =  1  eigenvalues  and  eigenvectors 
A  dimension  statement  in  the  main  program  for  sufficient  storage  space  of 
the  vectors  L,  M,   (utilized  by  EIGEN)  and  R   (utilized  by  MINV)  is  also 
necessary.   All  of  the  required  information  to  construct  the  table  of 
sensitivity  data  is  now  available  in  the  main  program.   Insertion  of  the 
proper  equations  and  output  statements  completes  the  modification  for 
sensitivity  analysis. 

A  listing  of  MINV  and  EIGEN  and  an  example  output  format  follows. 
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SUL3R0UT  IN  li    MINV 

purpose 

invert  a  katrix 

U  S  A  G  E 

CALL  M  I  N  V  (  A  ,  N  ,  D  ,  L  ,  M  ) 

DESCRIPTION  OF  P  A  :<  AM  E '1  rR  S  .   ,., 

A  -  INPUT  MATRIX,  DESTROYED  IN  COMPUTATION  AND  rtEj?L;iC* 

FESULTANT  INVERSE. 
N  -  ORDER  0?  MATRIX  A 
D  -  RESULTANT  DETERMINANT 
L  -  fcORK  VECTOR  OF  LEN.GTE  N 
H  -  WORK  VECTOR  CE  LENGTH  K 

REMARKS 

'•MATRIX  A  MUST  EE  A  GENERAL  MATRIX 

SUBROUTINES  AND  FUNCTION  SUEPROGRAMS  REQUIRED 
NONE 


METHOD 

TilE  STANDARD  GAUSS- JORDAN  METHOD  IS  USED.  THE  DETERMINE  .' 
IS  ALSO  CALCULATED.  A  DETERMINANT  OF  ZERO  INDICATES  IMT* 
THE  MATRIX  IS  SINGULAR. 


.  ■ 
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fcUU  ECU  I I  N  IS  Ml  N V ( A, N , D, L  ,  M) 
DIMENSION  Ml)  #L  O)  iB(1) 


15  A  DOUliLi?  PRECISION  VERSION  Of  THIS  ROUTINE  IS  DE:  "  F  I ,  'I 
C  IN  COLUMN  \    SHOULD  Br  REMOVED  FECM  THE  JCU/L;  PRECIS!  . 
STATEMENT  WHICH  FOLLOWS. 


DOUBLE  FFECISICN  A, D , EIG A / HOLD 


THE  C  MUST  ALSO  BE  REMOVED  FROM  DOUBLE  PRECISION  STATIMEMS 
APPEARING  IN  OTHER  ROUTINES  USEE  IN  CONJUNCTION  hTTfi  THIS 
ROUTINE. 


V 


THE  DOUBLE  PRECISION  VERSION  CE  THIS  SUBROUTINE  MUST  ALSO 
CONTAIN  DOUBLE  PRECISION  rORTEAN  FUNCTIONS.   ABS  IN  STAIIME1 
10  MUST  3z    CHANGED  TO  EABS. 


SEARCH  5CE  LARGEST  ELEMENT 


1=1.  G 
IK=-N 

EO  EC  K=1,N 

RK=KK+N 

L(N)=K 

I  (  N  )  =  K 

K  K  =  N  K  ♦  K 

3IGA=A  (KK) 

DC  20  J=K,N 

IZ=N*  (J-1) 

DC  20  I=K,S 

IJ=IZ+I 

Ir(  ABS(BIGi)- 

EIGA=A  (10) 

I(K)=I 

M  (K)  =J 

CONTINUE 


ABS(A(IJ)))  15,20,20 


INTERCHANGE  ROWS 

I=L(K) 

IP(J-K)  35,35,25 
IlrK-N 

DC  30  1=1, N 
KI=KI+N 
HOLD  =  -A  (KI) 
JI=KI-K+J 
A  (KI)=A  (01) 
A(JI)  =HCIE 

INTERCHANGE  COLUMNS 


I=K  (K) 

IE(I-K)  45,45,38 

JP=N*  (1-1) 

DO  40  J=1,N 

JK=NK*J 
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J  I  J  L>  ♦  J 
|OLD=-A  (JK) 
A  (JK)  =  A  (JI) 
A  (01)  =ii CLE 

divide  column  3y  minus  eivci  (value  of  pivot  ei 

CONTAINED    IN    i-IGA) 

tF(3IGA)     48,46,48 

D=  0 .  0 

ETUBN 

CO    5  5    1=1, N 

I?  (I-K)     50,55,50 

IK=NK+I 

A  (IK)  =  A  (IK)/  (-BIGA) 

CONTINUE 

SEDUCE    KA1SIX 

DC    65    1=1  ,N 

IK=NK+I 

|CLD=A  (IK) 

IJ=I-N 

DC     6  5    J=1,N 

IJ=IJ+N 

IF  (I-K)     6C,65,6C 

IE(J-K)     62,65,62 

KJ=IJ-I+K 

A  (IJ)  =HCLD*A(KJ)  +A(IJ) 

CONTINUE 

EIVIEI    BOft    BY    PIVOT 

KJ=K-K 

DO    7  5    J-1,.K 

IJ-=KJ  +  N 

IF(J-K)     7C,75,7C 
A (KJ)=A (KJ)/BIGA 

CONTINUE 

PRODUCT  OF  PIVOTS 
D=D*EIGA 

REPLACE  FIVCT  BY  5ECIPROCAL 

A  (KK)  =  1.0/EIGA 

CONTINUE 

FINAL  SO w  AND  COLUMN  INTERCHANGE 


MI 

.  r 

MI 

;•  t 


K  =  N 

K=  (K-l) 

IF(K)  150,150,105 

|=L(K) 

IF  (I-K)  12C,120,108 
|C=N*  (K-1  ) 

Ir=n*  (i-i) 

DO  110  J  =  1  ,N 
JK=JQ+J 


A-10 


I0LC  =  A  (JK) 

i  i = j  Li + j 

IvlK)=-ft  (J  I) 
l(JI)     =iiCLD 

[=:•;  (;<) 

f(J-K)      1CC,100,12£ 

:0     13  0    1  =  1,  N 

blD=A  (KI) 
I=KI-K+J 
i(KI)=-A  (JI) 
(JI)     =  HCIE 

I    10    ICC 
ETUBN 


.    r 


A-ll 


SUBROUTINE     EIGEK 


pukpo 


COMPUTE  EIGENVALUES  AND  EIGENVECTORS  C?  A  .; 
MATRIX 


usag: 

caii  eiges  (a,r,n,mv) 

description  cf  parameters 

A  -  ORIGINAL  MATRIX  (SYMMETRIC),  DESTROYED  IN  COMPUTATION 
RESULTANT  EIGENVALUES  ARE  DEVELOPED  IN  DIAGONAL  OF 


MATRIX  A  IN  DESCENDING  CrIEE. 
E  -  FESUL1ANT  MATRIX  Or  EIGENVECTORS  (STORED  COLUMNWISE, 

IN  SAME  SEQUENCE  AS  EIGESVALUES) 
N  -  ORDER  CF  MATRICES  A  AND  ?. 
MV-  INPUT  CODE 

0  COMPUTE  EIGENVALUES  AND  EIGENVECTORS 

1  COMPUTE  EIGENVALUES  ONLY  (R  NERD  NOT  DE 
EIKESSICNEE  EUT  MUSI  STILL  APPEAR  IN  CALLIN 
SEQUENCE) 

BEHAEKS 

ORIGINAL  MATRIX  A  MUST  BE  EEAL  SYMMETRIC  (STORAGE  MODS*1) 

MATRIX  A  CANNOT  EE  IN  TEE  2b?.c     ICCATICN  AS  MATRIX  F 

SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED 
NONE 

METHOD 

DIAGCNALI2ATICN  METHOD  ORIGINATED  EY  JACOEI  AND  ADAPTED 
BY  VCN  NEUMANN  FOR  LARGE  COMPUTERS  AS  FCUKE  IN  'MATHE MATICAi 
METHODS  FOR  DIGITAL  COMPUTERS',  EDITED  EY  A.  RALSTON  AND 
K.S.  KILE,  JOHN  KIL2Y  AND  SONS,  NEW  YORK,  '1962,  CHAPTER  7 
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SU3 ECU TINE  iilGEN  (A,R,N,MV) 
DIMENSION  A (1)  ,B (1) 


IF  A  DOUBLE  PRECISION  VESSICN  OF  THIS  ROUTINE  IS  DES 
C  IN  COLUMN  1  S  H  C  U  L  £  BE  E  E  K  C  V I E  5  5  C  M  1 K  E  D  C  0  3L I  P  R 
STATEMENT  51HICH  FOLLOWS. 

DOUBLE  PRECISION  A, R, ANORH, ANRKX ,THR, X, Y,SINX  ,SIN X2  ,CCS X , 
1  CCSX2,SINCS, RANGE 

Till  C  MUST  ALSO  BE  REMOVED  FRCK  DOUBLE  PRECISION  STATEMEN 

APPEARING  IN  OTHEB  ROUTINES  USE!  IN  CONJUNCTION  WITH 

ROUTINE . 

THE  DGUELE  PRECISION  VERSION  OF  THIS  SUBROUTINE  KDST  MSC 
CONTAIN  DOUBLE  PRECISION  FOETEAK  FUNCTIONS.   SORT  IN  STAT. 
40,  68,  75,  AND  78  MUST  EE  CHANGED  TC  ESCRT.   A3S  IN  STATi 
62  MUST  BE  CHANGED  TC  DABS.  THE  CONSTANT  IN  STATEMENT  5  S^ 
BE  CHANGED  TC  1.0D-12. 


~  , 


ul: 


GENERATE  IDENTITY  KATSIX 

RANGE-1 .  CE-6 
pF(MV-l)  10,25,10 

IQ-=-)l 

DO  2G  J=1,K 

I  Q  a  I C  ♦  N 

DO  2G  1=1, N 
IJ-IC*I 

R(IJ)=C.C 

IF(I-J)  20,15,20 
R  (IJ)  =  1  .  C 
CONTINUE 

COMPUTE  INITIAL  ANL  FINAL  NCEMS  (ANOEM  AND  ANORMX) 

ANOPM=G.C 

DC  35  1=1, N 

DO  35  J=I,N 

IF  (I-J)  30,35,30 

lA=I*(J*J-J)/2 

&NORM=ANCRK*A  (IA) *A  (IA) 

CONTINUE 

IF  (A  NORM)   16  5,165,40 

ANORH=1  .4  14*SQRT  (ANORM) 

ANRMX-ANORF*RANGE/FLCAT(N) 

INITIALIZE  INDICATCRS  AND  CC^I-UIE  THRESHOLD,  IHE 

IND=C 
THR=ANCRM 
THR=THR/FLOAT (N) 
L=1 
M=L  +  1 


CCMEUTE  SIN  AND  CCS 
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hQ=  («  *M-M)  /2 
LQ  =  (L*L-I)/2 

L  K  =  I  +  K  c 

IF(    AES(A  (LK))-THR)     130,65,65 

IND=1 

LL=L+LC 

M  M  =  M  ♦  M  C 

X=0.5*  (A  (LL)-A  (MM)  ) 

Y  =  -A  (LM)/    SCKT  (A  (IK)  »A(LK)  *X*X) 

I F  ( X )     7  0,75,75 

Y=-Y 

felNX=Y/    SCr-T(2.0*{1.0+(    SCRT  (1.0-Y*Y)  )  )  ) 

EINX2=SINX*SINX 

|CSX  =  SCH1  (  1.G-SINX2) 

lOSX2-COSX*COSX 

sixes  -SINX*CCSX 

ROTATE  I  AND  M  COLUMNS 
ILC'=N*  (1-1) 

1kc-n*  (m-1) 

DO  125  1=1,  N 

10=  (I*I-I)/2 

Ir(I-L)  80,115,30 

IF(I-M)  85,115,90 

Itv.  =  I+MC 

GO  TO  95 

IM  =  K  +  IQ 

IF(I-L)  1CC,1C5,1G5 

IL  =  I  +  LC 

30  10  11 C 

IL=L+IQ 

K  =  A  (XL)  *COSX-A  (IH)  *SINX 

I  (IM)=A  (II)  *SIKX+A(IK)  *CCSX 

k  (II)=X 

|F(MV-1)  12  0,125,120 

ILH=ILC+I 

IMR=IMQ+I 

l-S  (113)  *CCSX-E  (IMP)  *SINX 

|<IMR)=R (IIR) *SINX+B  (1KB) *C CSX 

3  (i lb)  =  x 

IONTINUE 

X  =  2 .  0 * 7i  (LM)  *SINCS 

|=A  (LL)  *COSX2+A  (MM)  *SINX2-X 

l=A  (LL) *SIXX2+A  (MM) *COSX2+X 

I  (LM)  =  (A  (LL)-A  (MM)  )  *SINCS+A  (LM)  *  (CCSX  2-SINX  2) 

&  (LL)  =Y 

&  (MM)  =X 

TESTS  ECS  CCMELETION 

TEST  FCE  £S  =  LAST  COLUMN 
IF  (M-N)  135,140,135 
30  TO  6  0 


TEST  FOE  L  =  SECCNE  FECK  LAST  CCIUKN 
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r;- u- (n- i) )    u5,iso,m5 

L = L  ♦  1 

GO   TC    5f- 

IF(INC-l)      160,  155,  160 

INC=C 

GC    TC    50 

CCMEARE  THRESHOLD  WITH  PINAL  NORK 

IF  (TER-ANREX)  165,165,45 

SORT  EIGENVALUES  ANE  EIGEKVEC1CES 

DO  185  I=1,N 

IQ=IC+N 

LL=I  +  (I*I-I)/2 

IC=N*  (1-2) 

DO  185  J-I,N 

JC=OC+N 

HX=J*  (J*J-J)/2 

IF  (A  (LL)-A  £KK))  170,185,185 

k=A(LL) 

A(LL)=A  (MK) 

A  (HM)  =  X 

IF(KV-l)  175,185,175 

DO  180  K=1,N 

ILR=I0+K 

Imr=jc«-k 

X=K  (ILR) 
l(ILK)=R  (1KB) 
I  (1MB)  =X 
CONTINUE 

RETURN 

END 
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TYPICAL  OUTPUT  FORMAT 
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ThZ  INVSESS  HZ5SIAS  na^IX  IS 

3.50  5'4  63  2u 

1.G95  VJOOO  . 


1.  0051*4000 

2.  3C.3S-4&Q  ) 


»-  c  "*■  »  t.      v.  >  • ;    .•  1  y 


32  6.  2'>.;S  -ill;":.  .C.MJ2 

■115.5432  2G9.:S69>7 


.00  00  0 


VECTCF    CF    OPTICAL    X 
1.00007  1.00017 


fclGNTIE    V  A:\IATiON    D2L     *    =  C.I  COCO 

LCWABL?    CH&irfs    IN    X     WKS    I  >i  HEP  rriEr  STL  Y    i-'OR    ASSI^; 

♦  soa  -   c   t  x  (i) 

C  '•    -  0  .  C 1  •". r-  a 


EI    V  it  Til  A  T I G 


|0«AaL2    CHAK3H3    IN    X    VAHISD    SlrtULl ABEOUSLY    FOP    ASSIGNEE    VAEIAT 


EI2VA 


1.39868 


).31772         0.63301 
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